Fit stomach content models

Author

Max Lindmark

Published

September 25, 2023

Load libraries

# Load libraries, install if needed
pkgs <- c("tidyverse", "tidylog", "devtools", "sdmTMB", "sdmTMBextra", "terra", "mapplots",
          "viridis", "visreg", "modelr", "future", "kableExtra") 

if(length(setdiff(pkgs,rownames(installed.packages()))) > 0){

    install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
  
  }

invisible(lapply(pkgs, library, character.only = T))

# Import some plotting functions
# Source code for map plots
# You need: # devtools::install_github("seananderson/ggsidekick") # not on CRAN; library(ggsidekick)
devtools::source_url("https://raw.githubusercontent.com/maxlindmark/pred-prey-overlap/main/R/functions/map-plot.R")
options(ggplot2.continuous.colour = "viridis")

# devtools::install_github("seananderson/ggsidekick") # not on CRAN 
library(ggsidekick)
theme_set(theme_sleek())

# Set path
home <- here::here()

# For crossvalidation: paralell processing
plan(multisession)

set.seed(99) 
# To load entire cache in interactive r session, do: 
# qwraps2::lazyload_cache_dir(path = paste0(home, "/R/main-analysis/03-fit-diet-models_cache/html"))

Read old and new data, combine!

df <- read_csv(paste0(home, "/data/clean/stomachs.csv")) |> 
  mutate(depth_sc = (depth - mean(depth))/sd(depth),
         year_f = as.factor(year),
         month_f = as.factor(month),
         ices_rect = as.factor(ices_rect),
         pred_length_sc = (pred_length - mean(pred_length)) / sd(pred_length)) |> 
  filter(pred_length > 10) |> 
  filter(FR_spr < 0.4) |> # Check this by year... 
  filter(FR_sad < 0.4) |> 
  filter(FR_her < 0.4)

# Read the prediction grid...
pred_grid <- bind_rows(read_csv(paste0(home, "/data/clean/pred_grid_(1_2).csv")),
                       read_csv(paste0(home, "/data/clean/pred_grid_(2_2).csv")))

# Scale with respect to data!
df |> group_by(month) |> summarise(n = n()) |> arrange(desc(n))
# A tibble: 8 × 2
  month     n
  <dbl> <int>
1     3  4620
2    11  1842
3    12  1508
4     2   839
5     4   332
6    10    34
7     6    32
8     5     2
pred_grid <- pred_grid |> 
  drop_na(oxy, temp, sal, depth) |> 
  filter(quarter == 4) |> # Not needed in theory for saduria...
  mutate(depth_sc = (depth - mean(df$depth)) / sd(df$depth),
         year_f = as.factor(year),
         month_f = as.factor(3),
         year = as.integer(year),
         pred_length_sc = 0,
         ices_rect = as.factor(ices_rect)) |> 
  droplevels()

Plot data!

df |> 
  pivot_longer(c(FR_spr, FR_her, FR_sad)) |> 
  ggplot(aes(year, value)) +
  geom_jitter(height = 0, alpha = 0.5) + 
  coord_cartesian(ylim = c(0, 0.1)) +
  facet_wrap(~name, ncol = 1) +
  geom_smooth()

Set up a mesh

mesh <- make_mesh(df, c("X", "Y"), cutoff = 6)

ggplot() +
  inlabru::gg(mesh$mesh) +
  coord_fixed() +
  geom_point(aes(X, Y), data = df, alpha = 0.2, size = 0.5) +
  annotate("text", -Inf, Inf, label = paste("n knots = ", mesh$mesh$n), hjust = -0.3, vjust = 3) + 
  labs(x = "Easting (km)", y = "Northing (km)")

Spatial and temporal cross validation, and AIC, to compare spatial vs non spatial models

5-fold spatial cross validation & AIC to compare two models: spatial fields or ices rectangle as random effects

# Set up spatial clusters
clust <- kmeans(df[, c("X", "Y")], 5)$cluster

# Plot
df$clust_id <- clust

plot_map +
  geom_point(data = df, aes(X*1000, Y*1000, color = as.factor(clust_id))) +
  scale_color_viridis(discrete = TRUE) + 
  geom_sf(size = 0.1)
Warning: Removed 33 rows containing missing values (`geom_point()`).
Warning: colourbar guide needs continuous scales.

# Sprat
# ICES as a random effects
spr_space_1_cv <- sdmTMB_cv(
  FR_spr ~ 0 + year_f + depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect),
  data = df,
  mesh = mesh,
  spatial = "off",
  family = tweedie(link = "log"),
  fold_ids = clust,
  k_folds = length(unique(clust))
)
Running fits with `future.apply()`.
Set a parallel `future::plan()` to use parallel processing.
# Spatial random field
spr_space_2_cv <- sdmTMB_cv(
  FR_spr ~ 0 + year_f + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
  data = df,
  mesh = mesh,
  spatial = "on",
  family = tweedie(link = "log"),
  fold_ids = clust,
  k_folds = length(unique(clust))
)
Running fits with `future.apply()`.
Set a parallel `future::plan()` to use parallel processing.
# Herring
# ICES as a random effects
her_space_1_cv <- sdmTMB_cv(
  FR_her ~ 0 + year_f + depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect),
  data = df,
  mesh = mesh,
  spatial = "off",
  family = tweedie(link = "log"),
  fold_ids = clust,
  k_folds = length(unique(clust))
)
Running fits with `future.apply()`.
Set a parallel `future::plan()` to use parallel processing.
# Spatial random field
her_space_2_cv <- sdmTMB_cv(
  FR_her ~ 0 + year_f + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
  data = df,
  mesh = mesh,
  spatial = "on",
  family = tweedie(link = "log"),
  fold_ids = clust,
  k_folds = length(unique(clust))
)
Running fits with `future.apply()`.
Set a parallel `future::plan()` to use parallel processing.
# Saduria
# ICES as a random effects
sad_space_1_cv <- sdmTMB_cv(
  FR_sad ~ 0 + year_f + depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect),
  data = df,
  mesh = mesh,
  spatial = "off",
  family = tweedie(link = "log"),
  fold_ids = clust,
  k_folds = length(unique(clust))
)
Running fits with `future.apply()`.
Set a parallel `future::plan()` to use parallel processing.
# Spatial random field
sad_space_2_cv <- sdmTMB_cv(
  FR_sad ~ 0 + year_f + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
  data = df,
  mesh = mesh,
  spatial = "on",
  family = tweedie(link = "log"),
  fold_ids = clust,
  k_folds = length(unique(clust))
)
Running fits with `future.apply()`.
Set a parallel `future::plan()` to use parallel processing.

10-fold temporal cross validation & AIC to compare two models: spatial fields or ices rectangle as random effects

# Set up temporal clusters

# This is not very straightforward! Since with temporal clusters I don't have much spatial overlap

# Note we remove year from the model because it can't be estimated for single years
# clust <- as.numeric(as.factor(df$year))
# 
# # Plot
# df$clust_id <- clust
# 
# df <- df |> mutate(clust_id = round(clust_id/2))
# 
# sort(unique(df$clust_id))
# 
# df <- df |> mutate(clust_id = ifelse(!year == 2021, 1, 2))
# 
# plot_map_fc +
#   geom_point(data = df, aes(X*1000, Y*1000, color = year)) +
#   scale_color_viridis() + 
#   facet_wrap(~clust_id) +
#   geom_sf(size = 0.1)
# 
# # Sprat
# # ICES as a random effects
# spr_temp_1_cv <- sdmTMB_cv(
#   FR_spr ~ depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect),
#   data = df,
#   mesh = mesh,
#   spatial = "off",
#   family = tweedie(link = "log"),
#   fold_ids = clust,
#   k_folds = length(unique(clust))
# )
# 
# # Spatial random field
# spr_temp_2_cv <- sdmTMB_cv(
#   FR_spr ~ depth_sc + s(pred_length_sc, k=3) + (1|month_f),
#   data = df,
#   mesh = mesh,
#   spatial = "on",
#   family = tweedie(link = "log"),
#   fold_ids = clust,
#   k_folds = length(unique(clust))
# )
# 
# 
# # Herring
# # ICES as a random effects
# her_temp_1_cv <- sdmTMB_cv(
#   FR_her ~ depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect),
#   data = df,
#   mesh = mesh,
#   spatial = "off",
#   family = tweedie(link = "log"),
#   fold_ids = clust,
#   k_folds = length(unique(clust))
# )
# 
# # Spatial random field
# her_temp_2_cv <- sdmTMB_cv(
#   FR_her ~ depth_sc + s(pred_length_sc, k=3) + (1|month_f),
#   data = df,
#   mesh = mesh,
#   spatial = "on",
#   family = tweedie(link = "log"),
#   fold_ids = clust,
#   k_folds = length(unique(clust))
# )
# 
# 
# # Saduria
# # ICES as a random effects
# sad_temp_1_cv <- sdmTMB_cv(
#   FR_sad ~ depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect),
#   data = df,
#   mesh = mesh,
#   spatial = "off",
#   family = tweedie(link = "log"),
#   fold_ids = clust,
#   k_folds = length(unique(clust))
# )
# 
# # Spatial random field
# sad_temp_2_cv <- sdmTMB_cv(
#   FR_sad ~ depth_sc + s(pred_length_sc, k=3) + (1|month_f),
#   data = df,
#   mesh = mesh,
#   spatial = "on",
#   family = tweedie(link = "log"),
#   fold_ids = clust,
#   k_folds = length(unique(clust))
# )

AIC

# ICES as a random effects
fit_spr_m1 <- sdmTMB(
  FR_spr ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect), 
  data = df,
  mesh = mesh,
  time = "year",
  time_varying = ~1,
  spatial = "off",
  spatiotemporal = "off",
  family = tweedie(link = "log"))
Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`
sanity(fit_spr_m1)
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
fit_her_m1 <- sdmTMB(
  FR_her ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect), 
  data = df,
  mesh = mesh,
  time = "year",
  time_varying = ~1,
  spatial = "off",
  spatiotemporal = "off",
  family = tweedie(link = "log"))
Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`
sanity(fit_her_m1)
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
fit_sad_m1 <- sdmTMB(
  FR_sad ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect), 
  data = df,
  mesh = mesh,
  time = "year",
  time_varying = ~1,
  spatial = "off",
  spatiotemporal = "off",
  family = tweedie(link = "log"))
Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`
sanity(fit_sad_m1)
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
# Spatial random field
fit_spr_m2 <- sdmTMB(
  FR_spr ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f), 
  data = df,
  mesh = mesh,
  time = "year",
  time_varying = ~1,
  spatial = "on",
  spatiotemporal = "off",
  family = tweedie(link = "log"))
Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`
sanity(fit_spr_m2)
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
fit_her_m2 <- sdmTMB(
  FR_her ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f), 
  data = df,
  mesh = mesh,
  time = "year",
  time_varying = ~1,
  spatial = "on",
  spatiotemporal = "off",
  family = tweedie(link = "log"))
Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`
sanity(fit_her_m2)
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✖ `ln_tau_G` standard error may be large
ℹ `ln_tau_G` is an internal parameter affecting `sigma_G`
ℹ `sigma_G` is the random intercept standard deviation
ℹ Try simplifying the model, adjusting the mesh, or adding priors
✖ `sigma_G` is smaller than 0.01
ℹ Consider omitting this part of the model
✔ Range parameter doesn't look unreasonably large
fit_sad_m2 <- sdmTMB(
  FR_sad ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f), 
  data = df,
  mesh = mesh,
  time = "year",
  time_varying = ~1,
  spatial = "on",
  spatiotemporal = "off",
  family = tweedie(link = "log"))
Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`

Check residuals and from selected models

# Summary
summary(fit_spr_m2)
Spatial model fit by ML ['sdmTMB']
Formula: FR_spr ~ 0 + depth_sc + s(pred_length_sc, k = 3) + (1 | month_f)
Mesh: mesh (isotropic covariance)
Time column: year
Data: df
Family: tweedie(link = 'log')
 
                coef.est coef.se
depth_sc            0.33    0.07
spred_length_sc     0.66    0.08

Smooth terms:
                    Std. Dev.
sds(pred_length_sc)     32.42

Random intercepts:
        Std. Dev.
month_f      0.65

Time-varying parameters:
                 coef.est coef.se
(Intercept)-1993    -6.19    0.42
(Intercept)-1994    -5.94    0.44
(Intercept)-1995    -7.15    0.40
(Intercept)-1996    -6.00    0.36
(Intercept)-1997    -7.00    0.43
(Intercept)-1998    -5.69    0.48
(Intercept)-1999    -6.08    0.35
(Intercept)-2000    -6.20    0.35
(Intercept)-2001    -5.35    0.37
(Intercept)-2002    -7.47    0.43
(Intercept)-2003    -4.62    0.37
(Intercept)-2004    -6.99    0.37
(Intercept)-2005    -5.83    0.35
(Intercept)-2006    -6.07    0.34
(Intercept)-2007    -6.23    0.34
(Intercept)-2008    -6.11    0.34
(Intercept)-2009    -5.77    0.35
(Intercept)-2012    -4.97    0.35
(Intercept)-2013    -5.57    0.34
(Intercept)-2014    -4.93    0.39
(Intercept)-2018    -6.22    0.42
(Intercept)-2021    -5.94    0.32

Dispersion parameter: 0.26
Tweedie p: 1.44
Matérn range: 19.81
Spatial SD: 1.20
ML criterion at convergence: -2157.438

See ?tidy.sdmTMB to extract these values as a data frame.
summary(fit_her_m2)
Spatial model fit by ML ['sdmTMB']
Formula: FR_her ~ 0 + depth_sc + s(pred_length_sc, k = 3) + (1 | month_f)
Mesh: mesh (isotropic covariance)
Time column: year
Data: df
Family: tweedie(link = 'log')
 
                coef.est coef.se
depth_sc           -0.13    0.09
spred_length_sc    -0.43    0.07

Smooth terms:
                    Std. Dev.
sds(pred_length_sc)      38.1

Random intercepts:
        Std. Dev.
month_f         0

Time-varying parameters:
                 coef.est coef.se
(Intercept)-1993    -8.97    0.50
(Intercept)-1994    -7.34    0.44
(Intercept)-1995    -6.75    0.36
(Intercept)-1996    -7.70    0.37
(Intercept)-1997    -5.99    0.34
(Intercept)-1998    -6.53    0.45
(Intercept)-1999    -6.92    0.31
(Intercept)-2000    -6.92    0.31
(Intercept)-2001    -6.82    0.33
(Intercept)-2002    -7.55    0.45
(Intercept)-2003    -5.51    0.29
(Intercept)-2004    -6.48    0.28
(Intercept)-2005    -7.45    0.30
(Intercept)-2006    -6.96    0.27
(Intercept)-2007    -7.22    0.28
(Intercept)-2008    -6.99    0.27
(Intercept)-2009    -6.39    0.30
(Intercept)-2012    -5.88    0.33
(Intercept)-2013    -5.74    0.28
(Intercept)-2014    -6.96    0.47
(Intercept)-2018    -7.66    0.42
(Intercept)-2021    -7.37    0.22

Dispersion parameter: 0.45
Tweedie p: 1.46
Matérn range: 17.38
Spatial SD: 1.41
ML criterion at convergence: -203.006

See ?tidy.sdmTMB to extract these values as a data frame.

**Possible issues detected! Check output of sanity().**
summary(fit_sad_m2)
Spatial model fit by ML ['sdmTMB']
Formula: FR_sad ~ 0 + depth_sc + s(pred_length_sc, k = 3) + (1 | month_f)
Mesh: mesh (isotropic covariance)
Time column: year
Data: df
Family: tweedie(link = 'log')
 
                coef.est coef.se
depth_sc           -0.44    0.15
spred_length_sc    -0.04    0.09

Smooth terms:
                    Std. Dev.
sds(pred_length_sc)     12.28

Random intercepts:
        Std. Dev.
month_f      0.64

Time-varying parameters:
                 coef.est coef.se
(Intercept)-1993    -8.64    0.93
(Intercept)-1994    -8.67    0.91
(Intercept)-1995    -8.37    0.87
(Intercept)-1996    -8.33    0.86
(Intercept)-1997    -8.05    0.86
(Intercept)-1998    -8.66    0.92
(Intercept)-1999    -8.34    0.84
(Intercept)-2000    -8.01    0.84
(Intercept)-2001    -7.30    0.84
(Intercept)-2002    -8.32    0.86
(Intercept)-2003    -7.84    0.84
(Intercept)-2004    -7.86    0.84
(Intercept)-2005    -7.38    0.83
(Intercept)-2006    -7.89    0.83
(Intercept)-2007    -8.54    0.84
(Intercept)-2008    -8.48    0.84
(Intercept)-2009    -8.37    0.85
(Intercept)-2012    -7.92    0.85
(Intercept)-2013    -8.24    0.84
(Intercept)-2014    -7.86    0.83
(Intercept)-2018    -7.44    0.79
(Intercept)-2021    -7.56    0.77

Dispersion parameter: 0.65
Tweedie p: 1.54
Matérn range: 83.70
Spatial SD: 2.59
ML criterion at convergence: -1144.977

See ?tidy.sdmTMB to extract these values as a data frame.
df |> filter(year == 2018) |> distinct(FR_sad)
filter: removed 8,944 rows (97%), 265 rows remaining
distinct: removed 250 rows (94%), 15 rows remaining
# A tibble: 15 × 1
    FR_sad
     <dbl>
 1 0.0368 
 2 0.0426 
 3 0.128  
 4 0      
 5 0.00219
 6 0.00151
 7 0.00400
 8 0.00939
 9 0.0218 
10 0.0194 
11 0.0740 
12 0.00336
13 0.00233
14 0.00254
15 0.0115 
# Residuals
spr_res <- mcmc_res <- residuals(fit_spr_m2, type = "mle-mcmc",
                                 mcmc_samples = sdmTMBextra::predict_mle_mcmc(fit_spr_m2,
                                                                              mcmc_iter = 201,
                                                                              mcmc_warmup = 200))

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.002663 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 26.63 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 28.5231 seconds (Warm-up)
Chain 1:                0.05815 seconds (Sampling)
Chain 1:                28.5812 seconds (Total)
Chain 1: 
df$spr_res <- as.vector(spr_res)

her_res <- mcmc_res <- residuals(fit_her_m2, type = "mle-mcmc",
                                 mcmc_samples = sdmTMBextra::predict_mle_mcmc(fit_her_m2,
                                                                              mcmc_iter = 201,
                                                                              mcmc_warmup = 200))

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.002514 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 25.14 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 308.387 seconds (Warm-up)
Chain 1:                1.77732 seconds (Sampling)
Chain 1:                310.165 seconds (Total)
Chain 1: 
Warning: There were 1 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
Warning: Examine the pairs() plot to diagnose sampling problems
df$her_res <- as.vector(her_res)

sad_res <- mcmc_res <- residuals(fit_sad_m2, type = "mle-mcmc",
                                 mcmc_samples = sdmTMBextra::predict_mle_mcmc(fit_sad_m2,
                                                                              mcmc_iter = 201,
                                                                              mcmc_warmup = 200))

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.002502 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 25.02 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 66.2415 seconds (Warm-up)
Chain 1:                0.112717 seconds (Sampling)
Chain 1:                66.3542 seconds (Total)
Chain 1: 
df$sad_res <- as.vector(sad_res)

# Plot all together
df |> 
  rename(Sprat = spr_res,
         Herring = her_res,
         Saduria = sad_res) |> 
  pivot_longer(c(Sprat, Herring, Saduria)) |> 
  ggplot(aes(sample = value)) +
  facet_wrap(~name) +
  stat_qq(size = 0.75, shape = 21, fill = NA) +
  stat_qq_line() +
  labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
  theme(aspect.ratio = 1)
rename: renamed 3 variables (Sprat, Herring, Saduria)
pivot_longer: reorganized (Sprat, Herring, Saduria) into (name, value) [was 9209x25, now 27627x24]

Plot conditional effects, random effects and spatial predictions

# Depth
vis_spr_dep <- visreg(fit_spr_m2, xvar = "depth_sc", plot = FALSE)
vis_her_dep <- visreg(fit_her_m2, xvar = "depth_sc", plot = FALSE)
vis_sad_dep <- visreg(fit_sad_m2, xvar = "depth_sc", plot = FALSE)

vis_dep <- bind_rows(vis_spr_dep$fit |> mutate(species = "Sprat"),
                     vis_her_dep$fit |> mutate(species = "Herring"), 
                     vis_sad_dep$fit |> mutate(species = "Saduria")) |> 
  mutate(var = "Depth (scaled)") |> 
  rename(x = depth_sc)

d_dep <- bind_rows(vis_spr_dep$res |> mutate(species = "Sprat"),
                   vis_her_dep$res |> mutate(species = "Herring"),
                   vis_sad_dep$res |> mutate(species = "Saduria")) |> 
  mutate(var = "Depth (scaled)") |> 
  rename(x = depth_sc)

# Month
vis_spr_mon <- visreg(fit_spr_m2, xvar = "month_f", plot = FALSE)
vis_her_mon <- visreg(fit_her_m2, xvar = "month_f", plot = FALSE)
vis_sad_mon <- visreg(fit_sad_m2, xvar = "month_f", plot = FALSE)

vis_mon <- bind_rows(vis_spr_mon$fit |> mutate(species = "Sprat"),
                     vis_her_mon$fit |> mutate(species = "Herring"), 
                     vis_sad_mon$fit |> mutate(species = "Saduria")) |> 
  mutate(var = "Month") |> 
  rename(x = month_f) |> 
  mutate(x = as.numeric(as.character(x)))

d_mon <- bind_rows(vis_spr_mon$res |> mutate(species = "Sprat"),
                   vis_her_mon$res |> mutate(species = "Herring"),
                   vis_sad_mon$res |> mutate(species = "Saduria")) |> 
  mutate(var = "Month") |> 
  rename(x = month_f) |> 
  mutate(x = as.numeric(as.character(x)))

# Predator length
vis_spr_len <- visreg(fit_spr_m2, xvar = "pred_length_sc", plot = FALSE)
vis_her_len <- visreg(fit_her_m2, xvar = "pred_length_sc", plot = FALSE)
vis_sad_len <- visreg(fit_sad_m2, xvar = "pred_length_sc", plot = FALSE)

vis_len <- bind_rows(vis_spr_len$fit |> mutate(species = "Sprat"),
                     vis_her_len$fit |> mutate(species = "Herring"),
                     vis_sad_len$fit |> mutate(species = "Saduria")) |> 
  mutate(var = "Predator length") |> 
  rename(x = pred_length_sc)

d_len <- bind_rows(vis_spr_len$res |> mutate(species = "Sprat"),
                   vis_her_len$res |> mutate(species = "Herring"),
                   vis_sad_len$res |> mutate(species = "Saduria")) |> 
  mutate(var = "Predator length") |> 
  rename(x = pred_length_sc)

vis <- bind_rows(vis_dep, vis_mon, vis_len)
vis_dat <- bind_rows(d_dep, d_mon, d_len)

# Plot!
ggplot(vis |> filter(!var == "Month"), aes(x = x, y = visregFit)) + 
  facet_grid(var ~ species) + 
  geom_point(data = vis_dat |> filter(!var == "Month"), aes(x = x, y = visregRes),
             alpha = 0.3, color = "gray50") +
  geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), alpha = 0.3, color = NA) +
  geom_line(color = "steelblue", linewidth = 1) + 
  labs(x = "Scaled variable", y = "Prediction") +
  NULL

ggsave(paste0(home, "/figures/conditional.pdf"), width = 17, height = 11, units = "cm")

# Now do month (categorial)
ggplot(vis |> filter(var == "Month"), aes(x = as.factor(x), y = visregFit)) + 
  facet_wrap(~ species) + # free grid?
  geom_jitter(data = vis_dat |> filter(var == "Month"), aes(x = as.factor(x), y = visregRes),
             alpha = 0.1, color = "gray70") +
  geom_errorbar(aes(ymin = visregLwr, ymax = visregUpr), alpha = 0.6, color = "steelblue",
                width = 0, linewidth = 1) +
  geom_point(color = "steelblue", size = 2) + 
  labs(x = "Month", y = "Prediction") +
  theme(aspect.ratio = 5/6) +
  NULL

ggsave(paste0(home, "/figures/conditional_month.pdf"), width = 17, height = 6, units = "cm")

Calculate indices

# Need to refit the models with the extra time argument. The reason we don't do that before is because it creates a mismatch in residual dimensions and data
# This is for interpolating between year using the random walk
extra_time <- pred_grid |> filter(!year %in% df$year) |> distinct(year) |> pull(year)
filter: removed 357,852 rows (76%), 113,862 rows remaining
distinct: removed 113,855 rows (>99%), 7 rows remaining
# Spatial random field
fit_spr_m2_extra <- sdmTMB(
  FR_spr ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f), 
  data = df,
  mesh = mesh,
  time = "year",
  time_varying = ~1,
  extra_time = extra_time,
  spatial = "on",
  spatiotemporal = "off",
  family = tweedie(link = "log"))

fit_her_m2_extra <- sdmTMB(
  FR_her ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f), 
  data = df,
  mesh = mesh,
  time = "year",
  time_varying = ~1,
  extra_time = extra_time,
  spatial = "on",
  spatiotemporal = "off",
  family = tweedie(link = "log"))

fit_sad_m2_extra <- sdmTMB(
  FR_sad ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f), 
  data = df,
  mesh = mesh,
  time = "year",
  time_varying = ~1,
  extra_time = extra_time,
  spatial = "on",
  spatiotemporal = "off",
  family = tweedie(link = "log"))

# Predict on grid, for indices and maps
pred_spr <- predict(fit_spr_m2_extra, newdata = pred_grid, return_tmb_object = TRUE)
pred_her <- predict(fit_her_m2_extra, newdata = pred_grid, return_tmb_object = TRUE)
pred_sad <- predict(fit_sad_m2_extra, newdata = pred_grid, return_tmb_object = TRUE)

# Make temporal index!
ncells <- filter(pred_grid, year == max(pred_grid$year)) |> nrow()
filter: removed 455,448 rows (97%), 16,266 rows remaining
index_spr <- get_index(pred_spr, area = rep(1/ncells, nrow(pred_spr$data)))
Bias correction is turned off.
It is recommended to turn this on for final inference.
index_her <- get_index(pred_her, area = rep(1/ncells, nrow(pred_her$data)))
Bias correction is turned off.
It is recommended to turn this on for final inference.
index_sad <- get_index(pred_sad, area = rep(1/ncells, nrow(pred_sad$data)))
Bias correction is turned off.
It is recommended to turn this on for final inference.
# Make long
index <- bind_rows(index_spr |> mutate(prey = "Sprat"),
                   index_her |> mutate(prey = "Herring"),
                   index_sad |> mutate(prey = "Saduria")) |> 
  mutate(Observed = ifelse(year %in% extra_time, "No", "Yes"))
mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'Observed' (character) with 2 unique values and 0% NA
# For comparison with data
df_sum <- df |>
  group_by(year) |>
  summarise(Sprat = mean(FR_spr),
            Herring = mean(FR_her),
            Saduria = mean(FR_sad)) |> 
  pivot_longer(c("Sprat", "Herring", "Saduria"), names_to = "prey", values_to = "mean_fr")
group_by: one grouping variable (year)
summarise: now 22 rows and 4 columns, ungrouped
pivot_longer: reorganized (Sprat, Herring, Saduria) into (prey, mean_fr) [was 22x4, now 66x3]
sort(unique(index$year))
 [1] 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007
[16] 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021
sort(unique(pred_grid$year))
 [1] 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007
[16] 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021
index |> 
  ggplot(aes(year, est, fill = Observed)) +
  geom_point(shape = 21, alpha = 0.7) +
  scale_fill_manual(values = c("white", "grey10")) +
  facet_wrap(~prey, ncol = 1, scales = "free") + 
  geom_errorbar(aes(ymin = lwr, ymax = upr), alpha = 0.4, width = 0) +
  geom_point(data = df_sum, aes(year, mean_fr, color = "Data"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
  labs(x = "Year", y = "Feeding ratio", color = "") +
  scale_x_continuous(breaks = seq_range(unique(df$year), by = 3)) +
  theme(legend.position = "bottom")

ggsave(paste0(home, "/figures/supp/index_ci.pdf"), width = 11, height = 20, units = "cm")

index |> 
  mutate(upr = ifelse(prey == "Sprat" & upr > 0.015, 0.015, upr)) |> 
  ggplot(aes(year, est)) +
  geom_point(alpha = 0.7) + 
  stat_smooth(method = "gam", formula = y~s(x, k=4), color = "steelblue") +
  facet_wrap(~prey, ncol = 3, scales = "free") + 
  labs(x = "Year", y = "Feeding ratio", color = "") +
  scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
  theme_sleek(base_size = 9) +
  theme(legend.position = c(0.9, 0.9))
mutate: changed 10 values (11%) of 'upr' (0 new NA)

ggsave(paste0(home, "/figures/index.pdf"), width = 17, height = 6, units = "cm")

Plot spatial predictions

spatial_preds <- bind_rows(pred_spr$data |> mutate(prey = "Sprat"),
                           pred_her$data |> mutate(prey = "Herring"),
                           pred_sad$data |> mutate(prey = "Saduria"))
mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA
# Plot spatial random effect
plot_map +
  geom_raster(data = spatial_preds |> filter(year == 2000),
              aes(X*1000, Y*1000, fill = omega_s)) +
  scale_fill_gradient2(name = "Spatial random field") +
  facet_wrap(~prey) +
  geom_sf(size = 0.1) + 
  theme_sleek(base_size = 9) + 
  theme(legend.position = c(0.05, 0.65)) +
  theme(legend.position = c(0.05, 0.7),
        legend.key.height = unit(0.6, "line"),
        legend.key.width = unit(0.2, "line"))
filter: removed 1,366,344 rows (97%), 48,798 rows remaining
Warning: Removed 2283 rows containing missing values (`geom_raster()`).

ggsave(paste0(home, "/figures/supp/diet_omega_s.pdf"), width = 17, height = 6, units = "cm")
Warning: Removed 2283 rows containing missing values (`geom_raster()`).
# Plot spatiotemporal predictions
plot_map_fc +
  geom_raster(data = spatial_preds |> filter(prey == "Sprat"), aes(X*1000, Y*1000, fill = exp(est))) +
  scale_fill_viridis(trans = "sqrt", name = "FR") +
  facet_wrap(~year) +
  geom_sf(size = 0.1)
filter: removed 943,428 rows (67%), 471,714 rows remaining
Warning: Removed 22069 rows containing missing values (`geom_raster()`).

ggsave(paste0(home, "/figures/supp/sprat_spatiotemporal_diet.pdf"), width = 17, height = 17, units = "cm")
Warning: Removed 22069 rows containing missing values (`geom_raster()`).
plot_map_fc +
  geom_raster(data = spatial_preds |> filter(prey == "Herring"), aes(X*1000, Y*1000, fill = exp(est))) +
  scale_fill_viridis(trans = "sqrt", name = "FR") +
  facet_wrap(~year) +
  geom_sf(size = 0.1)
filter: removed 943,428 rows (67%), 471,714 rows remaining
Warning: Removed 22069 rows containing missing values (`geom_raster()`).

ggsave(paste0(home, "/figures/supp/herring_spatiotemporal_diet.pdf"), width = 17, height = 17, units = "cm")
Warning: Removed 22069 rows containing missing values (`geom_raster()`).
plot_map_fc +
  geom_raster(data = spatial_preds |> filter(prey == "Saduria"), aes(X*1000, Y*1000, fill = exp(est))) +
  scale_fill_viridis(trans = "sqrt", name = "FR") +
  facet_wrap(~year) +
  geom_sf(size = 0.1)
filter: removed 943,428 rows (67%), 471,714 rows remaining
Warning: Removed 22069 rows containing missing values (`geom_raster()`).

ggsave(paste0(home, "/figures/supp/saduria_spatiotemporal_diet.pdf"), width = 17, height = 17, units = "cm")
Warning: Removed 22069 rows containing missing values (`geom_raster()`).
# Plot spatiotemporal predictions for a year and all species
plot_map +
  geom_raster(data = spatial_preds|> filter(year == "2000"), aes(X*1000, Y*1000, fill = exp(est))) +
  scale_fill_viridis(trans = "sqrt", name = "FR") +
  facet_wrap(~prey) +
  geom_sf(size = 0.1) + 
  theme_sleek(base_size = 9) + 
  theme(legend.position = c(0.05, 0.7),
        legend.key.height = unit(0.6, "line"),
        legend.key.width = unit(0.2, "line"))
filter: removed 1,366,344 rows (97%), 48,798 rows remaining
Warning: Removed 2283 rows containing missing values (`geom_raster()`).

ggsave(paste0(home, "/figures/supp/spatial_prey_prediction_2000.pdf"), width = 17, height = 6, units = "cm")
Warning: Removed 2283 rows containing missing values (`geom_raster()`).

Plot correlations between spatial overlap and population-level FR

# Join feeding ratio indicies (per capita) and overlap
cod_spr <- read_csv(paste0(home, "/output/cod_pel_sum_ovrlap.csv")) |>
  dplyr::select(year, cod_spr_ovr_tot) |> 
  rename(overlap = cod_spr_ovr_tot) |> 
  mutate(prey = "Sprat")
Rows: 27 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (5): year, cod_spr_ovr_tot, cod_spr_ovr_sd_tot, cod_her_ovr_tot, cod_her...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed one variable (overlap)

mutate: new variable 'prey' (character) with one unique value and 0% NA
cod_her <- read_csv(paste0(home, "/output/cod_pel_sum_ovrlap.csv")) |>
  dplyr::select(year, cod_her_ovr_tot) |> 
  rename(overlap = cod_her_ovr_tot) |> 
  mutate(prey = "Herring")
Rows: 27 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (5): year, cod_spr_ovr_tot, cod_spr_ovr_sd_tot, cod_her_ovr_tot, cod_her...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed one variable (overlap)

mutate: new variable 'prey' (character) with one unique value and 0% NA
cod_ben <- read_csv(paste0(home, "/output/cod_ben_sum_ovrlap.csv")) |>
  dplyr::select(year, cod_sad_ovr_tot) |> 
  rename(overlap = cod_sad_ovr_tot) |> 
  mutate(prey = "Saduria")
Rows: 54 Columns: 3
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (3): year, quarter, cod_sad_ovr_tot

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed one variable (overlap)

mutate: new variable 'prey' (character) with one unique value and 0% NA
overlap <- bind_rows(cod_spr, cod_her, cod_ben)

index_ovr <- index |>
  left_join(overlap, by = c("year", "prey")) |> 
  drop_na()
left_join: added one column (overlap)
           > rows only in x     6
           > rows only in y  (  0)
           > matched rows     108    (includes duplicates)
           >                 =====
           > rows total       114
drop_na: removed 6 rows (5%), 108 rows remaining
# Join in cod biomass data to calculate snapshot predation
cod_pred <- read_csv(paste0(home, "/output/pred_cod.csv")) |> 
  filter(quarter == 4)
Rows: 943428 Columns: 38
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (3): substrate, month_year, ices_rect
dbl (35): X, Y, year, lon, lat, depth, quarter, month, oxy, temp, sal, sub_d...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
filter: removed 471,714 rows (50%), 471,714 rows remaining
# Left join cod biomass density onto the grid of diet predictions
spatial_preds_dens <- left_join(spatial_preds, cod_pred |>
                                  dplyr::select(est, X, Y, year) |> 
                                  rename(est_codbiom = est))
rename: renamed one variable (est_codbiom)
Joining with `by = join_by(X, Y, year)`left_join: added one column (est_codbiom)
           > rows only in x           0
           > rows only in y  (        0)
           > matched rows     1,415,142
           >                 ===========
           > rows total       1,415,142
# Multiply local feeding ration with biomass density of cod
spatial_preds_dens <- spatial_preds_dens |> 
  mutate(pred = exp(est) * est_codbiom)
mutate: new variable 'pred' (double) with 1,415,142 unique values and 0% NA
# Summarise across years
spatial_preds_dens_sum <- spatial_preds_dens |> 
  group_by(prey, year) |> 
  summarise(tot_pred = sum(pred))
group_by: 2 grouping variables (prey, year)
summarise: now 87 rows and 3 columns, one group variable remaining (prey)
# Left_join with index data
index_ovr <- index_ovr |>
  left_join(spatial_preds_dens_sum, by = c("year", "prey")) |> 
  drop_na()
left_join: added one column (tot_pred)
           > rows only in x     0
           > rows only in y  (  6)
           > matched rows     108
           >                 =====
           > rows total       108
drop_na: no rows removed
# Check some regressions...
summary(lm(tot_pred ~ overlap, data = index_ovr |> filter(prey == "Sprat")))
filter: removed 81 rows (75%), 27 rows remaining

Call:
lm(formula = tot_pred ~ overlap, data = filter(index_ovr, prey == 
    "Sprat"))

Residuals:
   Min     1Q Median     3Q    Max 
-21096  -8859  -3394   4497  50120 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)    -3876      13500  -0.287    0.776
overlap        59012      41360   1.427    0.166

Residual standard error: 15530 on 25 degrees of freedom
Multiple R-squared:  0.0753,    Adjusted R-squared:  0.03831 
F-statistic: 2.036 on 1 and 25 DF,  p-value: 0.166
summary(lm(tot_pred ~ overlap, data = index_ovr |> filter(prey == "Herring")))
filter: removed 81 rows (75%), 27 rows remaining

Call:
lm(formula = tot_pred ~ overlap, data = filter(index_ovr, prey == 
    "Herring"))

Residuals:
    Min      1Q  Median      3Q     Max 
-4959.6 -3374.5 -2305.8   468.6 13247.2 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)     3900       4986   0.782    0.441
overlap         3804       9792   0.388    0.701

Residual standard error: 5332 on 25 degrees of freedom
Multiple R-squared:  0.006, Adjusted R-squared:  -0.03376 
F-statistic: 0.1509 on 1 and 25 DF,  p-value: 0.701
summary(lm(tot_pred ~ overlap, data = index_ovr |> filter(prey == "Saduria")))
filter: removed 54 rows (50%), 54 rows remaining

Call:
lm(formula = tot_pred ~ overlap, data = filter(index_ovr, prey == 
    "Saduria"))

Residuals:
    Min      1Q  Median      3Q     Max 
-4230.7 -2127.8   -73.2  1827.6  7535.0 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)     1140       1388   0.821    0.415   
overlap        12247       3512   3.487    0.001 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2719 on 52 degrees of freedom
Multiple R-squared:  0.1896,    Adjusted R-squared:  0.174 
F-statistic: 12.16 on 1 and 52 DF,  p-value: 0.001001
summary(lm(est ~ overlap, data = index_ovr |> filter(prey == "Sprat")))
filter: removed 81 rows (75%), 27 rows remaining

Call:
lm(formula = est ~ overlap, data = filter(index_ovr, prey == 
    "Sprat"))

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0049434 -0.0014152 -0.0004706  0.0008664  0.0080607 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept) -0.0002763  0.0024015  -0.115   0.9093  
overlap      0.0145854  0.0073571   1.983   0.0585 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.002762 on 25 degrees of freedom
Multiple R-squared:  0.1359,    Adjusted R-squared:  0.1013 
F-statistic:  3.93 on 1 and 25 DF,  p-value: 0.05851
summary(lm(est ~ overlap, data = index_ovr |> filter(prey == "Herring")))
filter: removed 81 rows (75%), 27 rows remaining

Call:
lm(formula = est ~ overlap, data = filter(index_ovr, prey == 
    "Herring"))

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0008822 -0.0005791 -0.0003672  0.0003922  0.0025145 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.0005373  0.0008206   0.655    0.519
overlap     0.0012624  0.0016116   0.783    0.441

Residual standard error: 0.0008775 on 25 degrees of freedom
Multiple R-squared:  0.02396,   Adjusted R-squared:  -0.01509 
F-statistic: 0.6136 on 1 and 25 DF,  p-value: 0.4408
summary(lm(est ~ overlap, data = index_ovr |> filter(prey == "Saduria")))
filter: removed 54 rows (50%), 54 rows remaining

Call:
lm(formula = est ~ overlap, data = filter(index_ovr, prey == 
    "Saduria"))

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0007962 -0.0005396 -0.0002388  0.0002035  0.0019474 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.0014897  0.0003855   3.865  0.00031 ***
overlap     -0.0003224  0.0009749  -0.331  0.74222    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.0007548 on 52 degrees of freedom
Multiple R-squared:  0.002098,  Adjusted R-squared:  -0.01709 
F-statistic: 0.1093 on 1 and 52 DF,  p-value: 0.7422
# Calculate correlations between FR, predation and overlap
cor <- plyr::ddply(index_ovr, c("prey"),
                   summarise,
                   cor_fr_ovr = round(cor(overlap, est), 2),
                   cor_pred_ovr = round(cor(overlap, tot_pred), 2)) |> 
  pivot_longer(c("cor_fr_ovr", "cor_pred_ovr")) |> 
  mutate(name = ifelse(name == "cor_fr_ovr", "Feeding ratio", "Total predation"))
summarise: now one row and 2 columns, ungrouped
summarise: now one row and 2 columns, ungrouped
summarise: now one row and 2 columns, ungrouped
pivot_longer: reorganized (cor_fr_ovr, cor_pred_ovr) into (name, value) [was 3x3, now 6x3]
mutate: changed 6 values (100%) of 'name' (0 new NA)
# Plot!
library(ggh4x)<

index_ovr |> 
  rename('Feeding ratio' = est, 
         'Total predation' = tot_pred) |> 
  pivot_longer(c('Feeding ratio', 'Total predation')) |> 
  ggplot(aes(overlap, value)) +
  geom_point() +
  #facet_wrap(name ~ prey, scales = "free") +
  ggh4x::facet_grid2(name ~ prey, scales = "free", independent = "y") +
  geom_smooth(aes(overlap, value), inherit.aes = FALSE, color = "steelblue",
              formula = y~s(x, k=10), method = "gam") +
  geom_text(data = cor,
            aes(label = paste("r=", value, sep = "")), x = -Inf, y = Inf, hjust = -.1, vjust = 2.5,
            inherit.aes = FALSE, fontface = "italic", size = 2.5, color = "tomato3") +
  labs(y = "Population-level feeding ratio", x = "Spatial overlap", color = "Year") +
  theme_sleek(base_size = 9) +
  theme(legend.position = "bottom", aspect.ratio = 1) +
  NULL
rename: renamed 2 variables (Feeding ratio, Total predation)
pivot_longer: reorganized (Feeding ratio, Total predation) into (name, value) [was 108x10, now 216x10]
Warning in library(ggh4x) < ggplot(pivot_longer(rename(index_ovr, `Feeding
ratio` = est, : longer object length is not a multiple of shorter object length
 [1]  TRUE  TRUE FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[13] FALSE  TRUE FALSE FALSE FALSE  TRUE FALSE  TRUE FALSE FALSE FALSE FALSE
[25] FALSE FALSE FALSE  TRUE  TRUE FALSE  TRUE FALSE FALSE
# Summarise the differences here, make fake facet grid free axis

ggsave(paste0(home, "/figures/fr_overlap_cor.pdf"), width = 17, height = 11, units = "cm")

TODO:

  • there are new data in the db that are not in stefans data also for the older period. find a way to not get duplicates.
  • update as I get more data
  • explore the spatial effects more
  • make a clean, straight script for all species that produces the plots for the report